Skip to content


This page gives hints on how to perform a Bethe-Salpeter calculation of neutral excitation energies and dielectric function with the ABINIT package.


Many-Body Perturbation Theory (MBPT [Onida2002]) defines a rigorous framework for the description of excited-state properties based on the Green’s function formalism. Within MBPT, one can calculate charged excitations (i.e. electron addition and removal energies), using for example Hedin’s GW approximation [Hedin1965] for the electron self-energy. In the same framework, neutral excitations (experimentally related to optical absorption spectra) are also well described through the solution of the Bethe-Salpeter equation (BSE [Onida2002],[Albrecht1998]). At present, the BSE represents the most accurate approach for the ab initio study of neutral excitations in crystalline systems as it includes the attractive interaction between electrons and holes thus going beyond the random-phase approximation (RPA) employed in the GW approximation.

The BSE implementation in ABINIT is discussed in the review article [Giantomassi2011].

BSE calculations are activated by setting optdriver= 99. The number of valence and conduction states in the e-h basis set is controlled by the variables bs_loband and nband, respectively. By default, the code solves the BSE in the so-called Tamm-Dancoff approximation (TDA) in which the coupling between resonant and anti-resonant transitions is ignored. A full BSE calculation including the coupling term can be done by setting bs_coupling= 1 in the input file. The variable bs_coulomb_term specifies the treatment of the screened interaction in the BSE kernel. By default, the code reads the W matrix from the SCR file generated by the GW code (optdriver= 3). Alternatively, one can model the spatial dependency of W with the model dielectric function proposed by Cappellini [Cappellini1993]. This option requires the specification of bs_coulomb_term= 21 and mdf_epsinf. The contribution due to the exchange term can be optionally excluded by using bs_exchange_term= 0. This is equivalent to computing the macroscopic dielectric function without local field effects. The frequency mesh for the macroscopic dielectric function is specified by bs_freq_mesh while zcut defines the complex shift to avoid the divergences due to the presence of poles along the real axis (from a physical standpoint, this parameter is related to the electron lifetimes). In addition, a scissor operator of energy mbpt_sciss can be used to correct the initial band energies and mimic the opening of the KS gap introduced by the GW approximation.

Three different algorithms for the solution of the BSE are implemented:

  • (1) Direct diagonalization,
  • (2) Haydock recursive algorithm,
  • (3) Conjugate gradient eigensolver.

The bs_algorithm input variable allows the user to choose among them. The Haydock algorithm [Haydock1980],[Benedict1998] is the recommended approach for the computation of optical spectra. The variables bs_haydock_niter and bs_haydock_tol control the Haydock iterative method.

Unfortunately, the Haydock solver does not give direct access to the eigenstates of the Hamiltonian, hence it cannot be used for the study of the excitonic wavefunctions. The conjugate gradient (CG) method employs standard iterative techniques [Payne1992] to compute the lowest eigenstates of the BSE Hamiltonian. This solver is more memory demanding than the Haydock approach since the eigenstates must be stored in memory, but it gives direct access to the excitonic states.

The CG algorithm should be preferred over the direct diagonalization especially when the number of eigenstates is much smaller than the size of the BSE Hamiltonian.

Note, however, that CG has been implemented only for TDA calculations (Hermitian matrices).

The most important results of the calculation are saved in five different files. The BSR file stores the upper triangle of the BSE resonant block in Fortran binary format (BSC for the coupling matrix). The HAYDR_SAVE file contains the coefficients of the tridiagonal matrix and the three vectors employed in the iterative algorithm. This file can be used to restart the calculation if convergence has not been achieved (related input variables gethaydock and irdhaydock). Finally, the macroscopic dielectric function with excitonic effects is reported in the EXC_MDF file while RPA_NLF_MDF and GW_NLF_MDF contain the RPA spectrum without local field effects obtained with KS energies and the GW energies, respectively.

Different schemes with different degrees of accuracy and computational load are available in ABINIT [Gillet2016] in order to improve the convergence rate of BSE calculations. They are are available in ABINITv8 with the bs_interp_* input variables.

In the non-spin-polarized case, spin-singlet as well as spin-triplet excitations are computed. Spin-polarized case is also available.





Selected Input Files




  • The bse tutorial tutorial on the Bethe-Salpeter Equation (BSE) deals with the computation of the macroscopic dielectric function of Silicon within the Bethe-Salpeter equation.